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Superconductor 

Shin-ichi Kaneko, Ken Matsuba, Muhammad Hafiz, Kcigo Yamasaki, Erika Kakizaki, Nobuhiko Nishida, 
Hiroyuki Takeya^, Kazuto Hirata\ Takuto Kawakami^*, Takeshi Mizushima^, Kazushige Machida^ 

Department of Physics, Tokyo Institute of Technology, Tokyo 152-8551, Japan, 
^National Institute for Materials Science, Tsukuha 305-0047, Japan, 
'^Department of Physics, Okayama University, Okayama 700-8530, Japan 

Quantized bound states at a vortex core are discretized in YNi2B2C. By using scanning 
tunneling spectroscopy witli an unprecedented 0.1 nm spatial resolution, we find and identify 
the localized spectral structure, where in addition to the first main peak with a positive low 
energy, a second subpeak coming fi^om the fourfold symmetric gap structure is seen inside 
the energy gap. Those spectral features are understood by solving the Bogoliubov-de Gennes 
equation for a fully three-dimensional gap structure. A particle-hole asymmetric spectrum at 
the core site and quantum oscillation in the spectra are clearly observed. 
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Quasiparticles (QPs) bound in a vortex core, such as 
Major ana particles, are expected to play a fundamen- 
tal role in various physical situations,^' ^ ranging from 
the static and dynamical properties of vortex matter to 
quantum computation where it is required to exchange 
QPs bound in a core to encode a qubit through non- 
Abclian statistics.^' In order to isolate and manipulate 
those QPs, we need to identify such particles experi- 
mentally and theoretically as a first step. Namely, we 
need to find the quantized energy levels for QPs using 
spectroscopic methods, such as scanning tunneling mi- 
croscopy and spectroscopy (STM/STS). This task is stiU 
far from being realized although several STM/STS ex- 
periments have probed the broad peak centered at the 
approximately zero-bias energy E = associated with 
the vortex bound states,''^' ^ namely, the so-called Caroli- 
de Gennes-Matricon (CdGM) state, ^ but those experi- 
ments did not resolve each discretized energy level fully 
quantum-mechanically. The origin of the difficulty in ob- 
serving a discretized core bound state is due to the fact 
that (A) T should be sufficiently low to satisfy the con- 
dition T < T^/ep {Tc is the superconducting transition 
temperature and is the Fermi energy). (B) The spa- 
tial resolution of STM/STS should be sufficiently fine to 
see the isolated bound state. Those two conditions have 
not been fulfilled in the exiting STM/STS experiments. 

In our previous STS study of the vortex bound state 
in YNi2B2C at 0.46K,^ we have observed broad peaks 
with electron-hole asymmetry in intensity in the tunnel- 
ing spectrum around the center of the vortex core and 
interpreted them in terms of a quantum regime vortex 
core where the wave function of the hole- like bound state 
starts from zero at the center and that of the electron-like 
bound state does from the maximum. In that study, the 
spatial resolution was about 1 nm and the average spec- 
tra in the 2x2 nm^ of the vortex core were measured. 
Here, we have performed an STM/STS experiment at 
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lower T with an unprecedented spatial resolution (~0.1 
nm) , expecting to clarify the quantum nature of vortex 
core states and to possibly observe new features such 
as the quantum oscillations of QPs in YNi2B2C where 
rc=15.6 K, i/c2=8.8 T, the mean free path ^=32 nm, 
the coherent length ,^=6 nm, and the penetration depth 
A=110 nm. We have also theoretically analyzed the ob- 
tained spectra to reveal that the STS data indeed exhibit 
crucial information on discretized quantum level struc- 
tures. 

A homemade scanning tunneling microscope can be 
operated at temperature down to 180 mK and at mag- 
netic fields of up to 6 T. We used a Pt-Ir alloy wire with 
a diameter of 0.3 mm that was mechanically sharpened 
in this experiment. Before cooling down the sample, we 
checked the tip condition over the Au film. STS mea- 
surement was conducted simultaneously with STM mea- 
surement. Using STM, we scanned the sample for a wide 
surface where the square lattice can be seen. We chose a 
suitable spot and size range for the STS measurement. 
Bias voltage and feedback current were set before the 
measurement. We measured the I-V characteristics of 
the sample and numerically calculated the differentiation 
of the I-V characteristics to find the differential conduc- 
tance. We normalized dl/dV with the slope of the I-V 
characteristics by fitting it linearly. Then the vortex was 
imaged by plotting the value of dl/dV at various biases. 

In YNi2B2C, a variety of experiments^"^ have shown 
on the existence of point nodes, such as angle-resolved 
specific heat^ and thermal conductivity.^ From an STS 
experiment, the tunneling spectrum was found to have 
an dependence near i? = in the superconducting 
gap at zero magnetic field. In order to explain the E^ 
dependence, the existence of point nodes with the asym- 
metric recovery and the need for 3D calculation were 
proposed.^'' Angle-resolved photoemission measurement 
shows a strong gap magnitude variation on the cylindri- 
cal Fermi surface in the 17th band at the X point of the 
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tetragonal Brillouin zone of this material.^ It is known 
that a two-dimensional c?-wave nodal gap does not yield 
an isolated vortex bound state inside the bulk gap be- 
cause all the low-energy wave functions are extended, 
leaking out from the nodal directions/^ in contrast to 
the isotropic s-wave gap.^^ Thus, it is not clear whether 
the point node case exhibits the CdGM like bound state. 

Our analytic strategy is to solve the Bogoliubov-de 
Gennes (BdG) equation in three dimensions (3D) for var- 
ious possible 3D gap functions, including the point node 
gap function on a cylindrical Fermi surface. The BdG for- 
malism deals with the quantum limiting regime beyond 
the quasiclassical approximation^'^^^^ in the sense of the 
accessibility to discretized eigenstates and their quantum 
oscillation with the order of lattice spacing. 

In general, the BdG equation in spin singlet supercon- 
ducting states is given as 
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where the single particle Hamiltonian Hq ~ 5{ri — 

r.) (^— ^y{^\) — ^F^ and the QP wave function 

Uvij) = [u,y(r), Wiy(r)]^ with the Fermi energy ep- The 
pair potential A(ri,r2) is decomposed into A(i2) in 
the center-of-mass coordinate R = {pi + p.)/^, where 
Pi = Xi + Di, and the gap function ip{r) in the relative 
coordinate r = ri — r.. Here, we consider the single vor- 
tex state described as A(J2) = Aq tanh(|i2|/^) exp(i(^), 
where the coherence length ^ — eF/i^Aokp) and the az- 
imuthal angle (j) = arctan[(a;i -I- X2)/{yi + 2/2)]- In this 
model, the wave number becomes a well-defined quan- 
tum number because of the translation symmetry along 
the 2-axis. The BdG equation (1) is reduced to the k^- 
resolved two-dimensional form 
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= E^,k,u^,kAPi)^ (2) 
where Ho = S{pi-p2) (-^ + V{pi)-ejP{k,)) and 



"V'- =d^+dy. The two-dimensional form of the Fermi en- 
ergy e'jP{kz) = ^{kp — k1) reflects the fc^-cross section 
of the Fermi surface. The index z/ G Z denotes the i^-th 
excited state of the fc^-resolved BdG equation, eq. (2). 
The gap function (pk^{p) = (5:^ / dk^dkyifke^^'P is ob- 
tained from ifk in the fe-space by performing the fast 
Fourier transform. We can solve the BdG equation for 
the arbitrary anisotropic function ipk. Then, we obtain 
the QP local density of states (LDOS) 

J^{p,E) ^E^kAp,E)^j: \u,,kAp)\^6AE^E,^kJ, 

kz y^kz 

with the Lorentzian function 5ri{z) = r]'^ /[z^ + rf]. r] is 
adjusted to compare the experimental data of the STS 
measurements below. The procedure of the numerical di- 
agonalization of the BdG equation is explained in ref. 17. 

It is considered that the gap function of YNi2B2C 
is suppressed at four points, where nesting vectors ex- 




Fig. 1. (Color online) Schematic view of the gap function and 
definition of the angle (p drawn on the cylindrical Fermi surface 
used in the calculation. 
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We parametrize the gap function as 
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"^fc = 2 ■f"'" + ^niinik;,) + [1 - (<5min(fcz)] COS40fe} , (3) 



where the angle (pk 
minimum (/3min(fcz) = 



= arctan(fcj./fcj,) and in-plane gap 
ifQ[ta'nh{kz/aki)]"'. This gap func- 



tion is equivalent to the s+g-wave gap function in the 
kz ~ plane, and the parameters a, n, and tpo signify 
how rapidly the point nodes at the fe^ = and ky = 
plane recover the full gap with increasing kz toward the 
Brillouin zone boundary at kz = ki. We choose the cylin- 
drical Fermi surface as the first approximation, where 
^jPikz) is constant. This is because we focus on the 17th 
Fermi surface around the X point of YNi2B2C, which 
does not close in the Brillouin zone as obtained by the 
band calculation.^^ Figure 1 shows a three-dimensional 
view of ip^ along with the cylindrical Fermi surface. As 
we will see later, the combination of the Fermi surface 
shape and the gap function along the kz direction is cru- 
cial for reproducing the STS data. We carry out the cal- 
culation within the quantum limit regime kp^ — 10. In 
the case of the YNi2B2C, the direction of the Fermi ve- 
locity vp a.t gap nodcs^^ is different by tt/A from that of 
the wave vector kp, so that the low-energy QP spectrum 
rotates by n/A from the isotropic case vp \\ kp}^ Thus, 
the spectrum along the [110] ([100])-axis in our experi- 
ments corresponds to that along the = (7r/4) in our 
calculations. 

Figure 2 shows the experimental (a) and theoretical 
(b) spectral evolutions along the directions (a) [110] and 
(b) = 0. In the experimental data shown in Fig. 2(a), 
it is evident to see that (i) the main peak at the core 
is situated on the positive energy side of about 0.3 meV 
with an asymmetric shape, (ii) The 2nd peak at the core 
site exists at the positive bias 1.5 meV, marked with 
the arrow in Fig. 2(a). (iii) The spatial modulated be- 
havior appears on the evolution ridge. We note that the 
average of the present STS data with a 0.1 nm spatial res- 
olution in a 2x2 nm^ well reproduces the previous STS 
results of tunneling spectra by Nishimori et al.^ These 
features (i)-(iii) are supported by the numerical calcula- 
tion shown in Fig. 2(b). In particular, in Fig. 2(b), the 
feature (iii) is interpreted as the quantum oscillations of 
the peak intensities associated with several lower energy 
excitations. We also find that (iv) quantum oscillation 
depends on either the = or 7r/4 direction. The finite 
energy spectrum near ii^/Ao^O.3 in the = direction 
extend longer than that in the = 7r/4 direction. This 
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Fig. 2. (Color online) (a) Differential conductance at r=190mK 
and f/=0.3T as a function of the sample bias obtained by STS 
measurements for [110] direction. The energy resolution is about 
0.1 meV in this measurment. The spacing between spectra is 0.11 
nm. The total path length is 11 nm. (b) The spectral evolutions 
of the LDOS J\f{E, r) for antinodal direction. The spatial scale is 
denoted by long arrows with a length ■ The short down arrows 
in the back indicate the shoulder structure at the core. It is seen 
that several low energy peaks near E = oscillate and evolve in 
their intensities as a function of the position from the core. 



angular dependence is also shown in the experimental 
data as the difference between the evolution of the [110] 
and [100] directions. These consistencies with the the- 
ory denote that the experimental data in Fig. 2 shows 
the quantum limiting behaviors of the vortex core bound 
states whose discreteness is specified by the energy scale 
A^/e^P. It is also found that the spherical or ellipsoidal 
Fermi surface smoothens the quantum oscillation of the 
LDOS. Hence, the cylindrical Fermi surface is necessary 
for the quantum oscillation. 

It is emphasized that the 2nd peak that appears at 
£'/Ao'~0.3, shown in Fig. 2(b), is a consequence of the 
rotational symmetry due to the point nodes and is never 
seen in the LDOS of isotropic s-wave superconductors.^^ 
In order to understand the physical origin of the 2nd 
peak, we decompose the spectrum into the fc^-resolved 
spectra as shown in Fig. 3(a), where the fc^-resolved spec- 
tra at the core are shown. The strong peak at E/ Aq = 0.1 
is associated with the lowest bound state iy=l. This iy=l 
peak becomes stronger and shifts to the positive side as 
kz increases because the gap becomes wider. Simultane- 
ously, the 2nd peak at i?/Ao = 0.3 grows and is distinctly 
separated from the main ^ = 1 peak when kz ^ 0. Fig- 
ure 3(b) shows the fc^-resolved eigenenergy of the BdG 
equation eq. (2). At /c^^O, where the in-plane gap min- 
imum is ipminikz = 0) = 0, the spectrum is continuous, 
and the eigenstates with E/Aq lower than the gap min- 
imum (fjnin become discrete with increasing fmin- The 
v = 5 eigenenergy is marked in Fig. 3 with filled symbols, 
which consists of the 2nd peak in the LDOS. 

The physical origin is explained as follows. As (pmin — > 
1, which is the limit of a two-dimensional isotropic s- 
wave case,^^ the branches of the eigenstates v approach 
the eigenstates quantized by the angular momentum 




Fig. 3. (Color online) (a) The fc^-resolved spectra A/j^^. (r = 0, -E) 
at the core. The solid (broken) line indicates inside (outside) of 
the in-plane gap minimum Aoipmin- It is seen that as increases 
the 2nd peak structure at E/Aq = 0.3 becomes split from the 
main park at E = 0.1. (b) The fc^-resolved eigenenergy E^ 
(symbols) and the gap minimum Aoi/Pmin (line). Filled symbols in 
the both panels denote the eigenenergy E,^ j. , the 5th excited 
state of the fcz resolved BdG equation (2). 



e = 1/2, 3/2, 5/2,- •• with the energy E'^+i.^.^^ = Iloq, 
where ujq ^ Aq / 2e'jP , and the wave function Ug^if.^{p) oc 
Jg_i{kFp) near the vortex core, where J,^{x) is the Bessel 
function^' The only £—1/2 state has a nonzero core- 
site wave function Ug_^i (p = 0) 7^ 0. When (pmin < 1, 
however, the fourfold rotationally symmetric pair poten- 
tial hybridizes with the eigenstates with the angular mo- 
mentum £ + 4n {neZ)}^ Thus, the ^=1/2 state can be 
hybridized to the £=l/2 + 4n states, which are realized 
as the v = An + l states described as 



(4) 



This 71' = contribution of the £ = 1/2 state in the wave 
function is the origin of the 2nd peak of the core- 
site LDOS. The peak arising from the i/ = 5 state is most 
enhanced because the £ = 9/2 state is the nearest-energy 
state from the ^=1/2 state. Although there are peaks 
from the i/ = — 3, 9 • ■ • states in the fcj-resolved core-site 
LDOS, these peaks cannot be seen in Fig. 2(a) or 2(b) 
with the kz integration. 

Next, we move to the two-dimensional profile of the 
LDOS. Figure 4 shows that the experimental and theo- 
retical two-dimensional LDOSs with different energies. It 
is known that the LDOS has a long tail along the nodal 
direction (/) = 7r/4 near E = 0; as the energy increases, the 
direction of the tails rotates to = in the case of the 
fourfold rotationally symmetric gap function at low en- 
ergies. The corresponding theoretical LDOS is shown in 
Figs. 4(c) and 4(d) where clear two-dimensional quantum 
oscillations are observed. Although in Figs. 4(a) and 4(b) 
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Fig. 4. (Color online) Experimental differential conductance 
dl/dV with the sample biases E = 0.0 meV (a) and E = 2.1 
meV (b) at T=180 mK and H=0.3 T. The corresponding two- 
dimensional LDOS map Af(E = 0.1Ao,p) (c) and J\f{E = 
0.36Ao,p). The panel (e) shows the schematic figure of the QP 
paths with the large contribution when the vortex core exists at 
the origin. 

those two-dimensional oscillations are hardly visible ow- 
ing to the lack of spatial resolution in this case, the quan- 
tum oscillations are well anticipated from Fig. 2(a) where 
the linear line cut of the two-dimensional LDOS with the 
finest spatial resolution is displayed. On the other hand, 
the depleted spectral weight on the vortex center and the 
surrounding four peaks observed in Fig. 4(b) are not re- 
produced in the theoretical LDOS in Fig. 4(d), because of 
the technical difficulty in the numerical diagonalization. 

The two-dimensional oscillatory behaviors are under- 
standable with the schematic picture of the QP path 
shown in Fig. 4(e). In the fourfold rotationally symmetric 
gap function, the low-energy QPs move from the (/) = 7r/4 
direction to the 0=^/4 direction around the vortex 
core.^'^"^^ The interference between the wave functions 
of the countercurrent QPs leads to a standing wave, that 
is, the oscillation of the LDOS. When the energy is E^O 
and the distance is ^ in Fig. 4(e), the QP "1" 
moving from A to B and the QP "2" moving from D 
to E interfere with each other. The resulting wavefront 
direction of the standing wave is = ±7r/4, as shown in 
Fig. 4(c). On the other hand, when the energy E becomes 



larger, the distance between paths becomes as large as 
lo^kp^ with the growth of the impact parameter. The 
large Iq allows the formation of the interference pattern 
by the QP "1" moving from B to C and the QP "3" mov- 
ing from F to B in Fig. 4(e). Thus, the LDOS oscillation 
pattern appears along the ^ = 0, 7r/2 direction in the 
high-energy region E/Aq = 0.36, as shown in Fig. 4(d). 
The change in the orientation of the LDOS is clearly ob- 
served in the experimental results in Fig. 4. The obser- 
vation of the rotation of the oscillation pattern remains 
a future problem. 

We have performed STM/STS experiments on 
YNi2B2C under an unprecedented spatial resolution of 
~0.1 nm at low temperatures. Through the full quantum 
mechanical analysis based on the Bogoliubov-de Gennes 
equation for three-dimensional gap structures, we have 
succeeded in identifying at the core position the quan- 
tized bound state at a positive energy as the first peak 
with the particle-hole asymmetric shape. We have also 
found the second peak induced by the fourfold symmetric 
gap structure. This is the first experimental observation 
of clearly discretized Caroli-de Gcnnes-Matricon vortex 
bound states. We conclude that, in reciprocal space, the 
point node positions correspond to the continuum spec- 
tra with a zero energy, while the gapped positions give 
rise to discrete spectra. Note that the energy of first peak 
at the vortex core center 6E ^ 0.3 meV in Fig. 2(a) 
is much larger than A'^/ep, where the Fermi energy is 
€f^^ eV. However, if we estimate that as SE= Ao/fci?^o, 
where the experimentally obtained value is fe^Co ~ 10, 
the resulting SE is reasonable. One of the reasons for the 
discrepancy between these estimates is that the realis- 
tic dispersion relations of YNi2B2C differ from the free 
electron model used in our analysis. In order to achieve 
quantitative estimates, one should take into account the 
realistic dispersion relations. 
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